home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / DPSTART / TEST1C / MYFIRSTM.C < prev   
Encoding:
C/C++ Source or Header  |  1996-04-18  |  2.0 KB  |  41 lines
  1. #include <iostream.h>       // required for input/output
  2. #include <Arrays_real.h>    // required for array functionality
  3. #include <CurvPlot.h>       // required for curve plotting
  4.  
  5. int main(int nargs, const char** args)
  6. {
  7.   initDIFFPACK(nargs, args); // not strictly necessary here, but a good habit
  8.   cout<<"Give number of solution points: ";
  9.   int n;            // declare an integer n reflecting the number of divisions
  10.   cin >> n;         // read n from standard input (cin)
  11.   real h=1.0/(n+1); // 1/(n+1) gives integer division (=0) so .0 is important!
  12.   MatTri(real) A(n);// Declaration of a tridiagonal matrix with n rows
  13.   Vec(real) b(n);   // Declaration of a vector of length n.
  14.   Vec(real) u(n);   // Declaration of a vector to hold the numerical solution.
  15.  
  16.   // --- Set up matrix A and vector b ---
  17.   A.fill(0.0);      // set all entries in A equal to 0.0
  18.   b.fill(0.0);      // set all entries in b equal to 0.0
  19.   for(int i=1;i<=n;i++)    // i++ means i=i+1
  20.     {
  21.       if(i != 1)   A(i,-1) = 1;   // init sub diagonal
  22.       if(i != n)   A(i,1)  = 1;   // init super diagonal
  23.                    A(i,0)  =-2;   // init main diagonal
  24.                    b(i)    =-1*h*h;
  25.       if (i==n)    b(i)    = b(i)-1;
  26.     }
  27.   A.printAscii(cout,"A before LU-fact.");       // print matrix to the screen
  28.   b.printAscii(cout,"right hand side");         // print vector to the screen
  29.   A.factLU();                                   // Make LU-factorisation
  30.   A.printAscii(cout,"A after LU-fact.");        // print matrix w/header
  31.   A.forwBack(b,u);                              // Solve Au=b
  32.   // in version 1.0: A.forwBackLU(b,u);         // Solve Au=b
  33.   cout<<"\n \n x           numerical          exact:\n"; // \n is newline
  34.   for(i = 1; i <= n; i++)
  35.     cout<<oform("%4.3f        %8.5f         %8.5f \n",i*h,u(i),i*h*(3-i*h)/2);
  36.   // this one requires the gnuplot program to be installed:
  37.   popUpPlot( h,  1-h,    u,    "gnuplot","700x500+0+0","solution","u");
  38.   //        xmin,xmax,solution, program , plot window , title    ,func.name
  39.   return 0;
  40. }
  41.